package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T32;
import scu.maqiang.fes.FES2T3PE;
import scu.maqiang.fes.FES2T62;
import scu.maqiang.fes.FES2T6PE;
import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class SOTS2DT36PE {
	public SOTS2DT36PE(FES2T3PE fsT3, FES2T6PE fsT6) {
		this.fsT3 = fsT3;
		this.fsT6 = fsT6;
	}

	public void elasticityCellSolution(MatrixFunc conFunc, double[][][] param, BVPType tp, double[][] FOCS,
			double[][] D, double[][] SOCS) {
		SRMatrix K = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(conFunc, param, tp, K);
		double[][] RHS = new double[5][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(conFunc, param, SOTS2DT36PE::elementPEirstOrder, tp, RHS);
		fsT3.applyBC_MBN(K, Direct.All, 1, 2, 3, 4);
		MVO.fill(FOCS, 0.0);
		NewIterSSolver solver = new NewIterSSolver(K);
		for(int i = 0; i < 5; i++) {
			fsT3.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4);
			solver.PCGSSOR(RHS[i], FOCS[i], 1.75, 1); // 11
		}

		computeHomogenizedPiezoelectricity(conFunc, param, FOCS, D);

		RHS = new double[6][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(conFunc, param, FOCS, D, SOTS2DT36PE::elementPESecondOrder, tp, RHS);
		MVO.fill(SOCS, 0.0);
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[2], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[3], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[4], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[5], Direct.All, 0.0, 1, 2, 3, 4);
		solver.PCGSSOR(RHS[0], SOCS[0], 1.75, 1);
		solver.PCGSSOR(RHS[1], SOCS[1], 1.75, 1);
		solver.PCGSSOR(RHS[2], SOCS[2], 1.75, 1);
		solver.PCGSSOR(RHS[3], SOCS[3], 1.75, 1);
		solver.PCGSSOR(RHS[4], SOCS[4], 1.75, 1);
		solver.PCGSSOR(RHS[5], SOCS[5], 1.75, 1);
	}

	public static void elementPEirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		double[][] B = new double[6][3];
		B[0][0] = B[1][2] = -grad[0][0] * area;
		B[2][0] = B[3][2] = -grad[1][0] * area;
		B[4][0] = B[5][2] = -grad[2][0] * area;
		B[1][1] = B[0][2] = -grad[0][1] * area;
		B[3][1] = B[2][2] = -grad[1][1] * area;
		B[5][1] = B[4][2] = -grad[2][1] * area;

		eleXY[0] = MVO.matmul(B, coef[0]);
		eleXY[1] = MVO.matmul(B, coef[1]);
		eleXY[2] = MVO.matmul(B, coef[2]);
	}

	public void computeHomogenizedPiezoelectricity(MatrixFunc conFunc, double[][][] param, double[][] FOCS, double[][] D) {
		double[][] grad = new double[3][2];
		double[][] xy = new double[3][2];
		Mesh2T3 mesh = (Mesh2T3) fsT3.getMesh();
		double[][] N = new double[3][6];
		double[][] B = new double[3][6];
		double wholeArea = 0.0;
		MVO.fill(D, 0.0);
		for (int k = 0; k < mesh.getNt(); k++) {
			int[] ele = mesh.getElement(k);
			mesh.getCoordPerElement(ele, xy);
			int[] dof = fsT3.getDof(ele);
			MVO.getSlice(FOCS[0], dof, N[0]);
			MVO.getSlice(FOCS[1], dof, N[1]);
			MVO.getSlice(FOCS[2], dof, N[2]);
			double area = CG.gradient2DT3(xy, grad);

			B[0][0] = B[2][1] = grad[0][0];
			B[0][2] = B[2][3] = grad[1][0];
			B[0][4] = B[2][5] = grad[2][0];
			B[1][1] = B[2][0] = grad[0][1];
			B[1][3] = B[2][2] = grad[1][1];
			B[1][5] = B[2][4] = grad[2][1];
			double[] ip = mesh.getinteriorPointInElement(xy);
			double[][] con = conFunc.action(ip, mesh.getElementLabel(k), param);
			for (int i = 0; i < 3; i++) {
				for (int j = 0; j < 3; j++) {
					D[i][j] += (con[i][j] + MVO.dot_product(con[i], MVO.matmul(B, N[j]))) * area;
				}
			}
			wholeArea += area;
		}
		MVO.mul(D, 1.0 / wholeArea);
	}

	public static void elementPESecondOrder(int[] ele, double[][] xy, double[][] con, double[][] FOCS, double[][] D,
			BVPType tp, double[][] eleXY) {
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		double area3 = area / 3.0;
		double[][] B = new double[3][6];
		B[0][0] = B[2][1] = grad[0][0];
		B[0][2] = B[2][3] = grad[1][0];
		B[0][4] = B[2][5] = grad[2][0];
		B[1][1] = B[2][0] = grad[0][1];
		B[1][3] = B[2][2] = grad[1][1];
		B[1][5] = B[2][4] = grad[2][1];
		double[][] BT = MVO.transpose(B);

		double[][] w = { { area3, 0.0 }, { 0.0, area3 }, { area3, 0.0 }, { 0.0, area3 }, { area3, 0.0 },
				{ 0.0, area3 } };
		double[][] wT = MVO.transpose(w);

		double[] C = new double[2];
		double[] part1 = null;
		double[][] con2 = new double[2][];
		double[][] con2T = null;
		double[] part2 = null;
		
		for(int k = 0; k < 3; k++) {
			C[0] = con[0][k] - D[0][k] + MVO.dot_product(con[0], MVO.matmul(B, FOCS[k]));
			C[1] = con[2][k] - D[2][k] + MVO.dot_product(con[2], MVO.matmul(B, FOCS[k]));
			part1 = MVO.matmul(w, C);
			
			con2[0] = con[0];
			con2[1] = con[2];
			con2T = MVO.transpose(con2);
			part2 = MVO.matmul(MVO.matmul(MVO.matmul(BT, con2T), wT), FOCS[k]);

			for (int i = 0; i < 6; i++) {
				eleXY[2 * k][i] = part1[i] - part2[i]; // 111
			}

			// 121
			C[0] = con[2][k] - D[2][k] + MVO.dot_product(con[2], MVO.matmul(B, FOCS[k]));
			C[1] = con[1][k] - D[1][k] + MVO.dot_product(con[1], MVO.matmul(B, FOCS[k]));
			part1 = MVO.matmul(w, C);

			con2[0] = con[2];
			con2[1] = con[1];
			con2T = MVO.transpose(con2);
			part2 = MVO.matmul(MVO.matmul(MVO.matmul(BT, con2T), wT), FOCS[k]);
			for (int i = 0; i < 6; i++) {
				eleXY[2 * k + 1][i] = part1[i] - part2[i]; // 121
			}
		}

	}

	public void assemblePESOTSSolution(FES2T32 fsFine, double[] T_0, double[] T_1, double[] T_2,
									   FES2T62 fsHomo, double[] T_Homo, FES2T32 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		
	}

	FES2T3PE fsT3;
	FES2T6PE fsT6;

}
